Sommerfeld Enhancement of DM Annihilation: 
Resonance Structure, Freeze-Out and CMB Spectral Bound 
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In the last few years there has been some interest in WIMP Dark Matter models featuring a 
velocity dependent cross section through the Sommerfeld enhancement mechanism, which is a non- 
relativist ic effect due to massive bosons in the dark sector. In the first part of this article, we find 
analytic expressions for the boost factor for three different model potentials, the Coulomb potential, 
the spherical well and the spherical cone well and compare with the numerical solution of the Yukawa 
potential. We find that the resonance pattern of all the potentials can be cast into the same universal 
form. 

In the second part of the article we perform a detailed computation of the Dark Matter relic density 
for models having Sommerfeld enhancement by solving the Boltzmann equation numerically. We 
calculate the expected distortions of the CMB blackbody spectrum from WIMP annihilations and 
compare these to the bounds set by FIRAS. We conclude that only a small part of the parameter 
space can be ruled out by the FIRAS observations. 



I. INTRODUCTION 

The most probable candidate for Dark Matter is ar- 
guably the WIMP, a Weakly Interacting Massive Par- 
ticle. This is due to what is often called the 'WIMP 
miracle'; new physics is required at the TeV scale when 
unitarity breaks down in the Standard Model, and some 
symmetry should protect the proton from decaying too 
fast through new diagrams involving these new degrees 
of freedom. The symmetry may ensure that at least one 
of the new particles is stable on cosmological timescales, 
while the annihilation cross section is typically within an 
order of magnitude of what is required to explain the 
observed relic abundance. 

In the last few years there have been some observa- 
tions p]-[9| pointing to anomalous emission within the 
galaxy. It is possible, that these observations may 
be explained by WIMPs annihilating in the galactic 
halo pjj-[l4|. Common to all these suggestions is that 
the necessary annihilation rate requires an annihilation 
cross section a few orders of magnitude larger than the 
one consistent with a thermal relic density, (cn;) TH c± 
3 • l(T 26 cm 3 /s. 

This apparent discrepancy, as well as the necessary 
suppression of hadronic final states, has been suggested 
to stem from new GeV-mass force carriers in the dark 
sector [l5|, [l6| which leads to Sommerfeld enhancement 
of the Dark Matter pair annihilation cross section. Al- 
though the anomalies mentioned above may turn out to 
be explained by normal astrophysical processes, there 
may still be force carriers in the dark sector which are 
light compared to the mass of the Dark Matter particle, 
so the mechanism of Sommerfeld enhancement is inter- 
esting and generic. Sommerfeld enhancement applied to 
WIMP annihilations in general, not related to the cosmic 



ray excess, was first studied in [lj 
and (H. 



hj and later in m 



II. SOMMERFELD ENHANCEMENT 

Sommerfeld enhancement is a consequence of having 
light force carriers mediate an attractive interaction in 
the dark sector, which creates a Yukawa-potential for the 
WIMPs. In this work we consider only s-wave annihila- 
tion. We write this potential as 
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where A is the coupling parameter, is the mass of 
the force carrier and r is the relative distance between 
the WIMPs. Since the potential is only a function of 
the relative distance, it only enters in the reduced one- 
particle Schrodinger equation for the relative motion. In 
spherical coordinates this is 
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where \i = m x /2 is the reduced mass. We are interested 
in the probability density of ipk when r = 0, since the 
boost factor for the s-wave case is given by Sk = 1^(0) | 2 , 
when tpk is normalized to the asymptotic form 
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This form of the boost factor can be derived [16| from a 
simple argument: Annihilations are assumed to proceed 
by a delta function interaction, so the rate of this process 
must be proportional to the norm squared of the reduced 
wave function at zero separation, a ex |^(r = 0)| 2 . De- 
noting the V = wave function by ifj° , we must have 
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When the potential is rotationally symmetric, we know 
that the solution of equation (j2j) becomes invariant under 
rotations around the axis of the incoming particle, and 
hence we can expand ip in products of Legendre poly- 
nomials and a radial wave function R^i. Only Rko is 
non-zero at the origin, so we define Xk = rRko, which 
leads to the equation 

1 _ ( ^^-m^r _ a2\^^ /gj 

(6) 
(7) 

where we introduced the dimensionless distance x = m^r 
and the ratio of the two masses / = m ( f ) /m x . Equation 
©, together with the boundary conditions 




x(o) = o 



X —> sin ( —x + 5 ] , as x — >> oo. 



(8a) 
(8b) 



defines the problem. 



A. Numerical solution 

In principle we should solve the boundary value prob- 
lem (|8j) numerically, but since the Schrodinger equation 
is linear, we can exchange the condition (|8b|) by setting 
the derivative of x to unity at x = and solve what is 
now an initial value problem for x- The solution to the 
original problem is just given by x — x/A where A is 
the asymptotic amplitude of x- The enhancement factor 
is then given by 



Sk = 
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A/3 
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In figure [T] we have plotted the Sommerfeld boost factor 
for a velocity of 150km/s for the whole (a, /) -parameter 
space. The vertical lines in figure [1] are ragged by thin, 
straight lines of resonance which have, as we shall see 
later, a rather large impact on the freeze out process. 

To calculate the amplitude A in ©, we must know 
when our waveform x nas reached its asymptotic form. 
This must happen when the Yukawa potential becomes 
much less than the kinetic energy of the particle, so we 
find the position of equality: 
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FIG. 1: Sommerfeld boost factor in log 10 for a relative velocity 
of 150km/s. 

where W is the Lambert W-function. To calculate the 
amplitude we use the asymptotic waveform (|8b|) : 



~ a ■ (P r 
X = Asm —x + o 

\f 

— = — A. cos \ -x + 6 
dx f \f 



(11a) 
(lib) 



A =\ X 2 + I ~s~r~ 



(3 dx 



which leads to this expression for the amplitude A: 

(12) 

We evolve our waveform x to 1.5x range , and after this 
point we calculate A at each succeeding point until A has 
converged. 



B. Coulomb potential 

In the limit of a massless force carrier, the Yukawa 
potential becomes a Coulomb potential. In this case, the 
enhancement factor can be calculated in closed form, 
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by solving the Schrodinger equation in terms of hyper 
geometric functions. The derivation of equation ([T3]) can 
be found in section [XI We want to stress, that the ap- 
plicability of this limit does not only depend on / <C 1, 
but is also dependent on a and j3. The correct question 
to ask is whether or not the scattering will take place in 
a regime in which the potential looks like a Coulomb po- 
tential. We want the range of the potential x range from 
([TO]) , to be less than something like ln(2). This would 
keep the exponential in the Yukawa potential at order 1 
during the interaction, so we should be safe in assum- 
ing a Coulomb potential instead. (However, this does 
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FIG. 2: Boost factor in log 10 for the Coulomb potential in 
log 10 for a relative velocity of 150km/s. 

not take into account the non-perturbative effects of res- 
onances.) The boost factor ([T3]) is plotted in figure [2j If 
we compare with the Sommerfeld-case in figure [TJ we can 
see that the lower left part of the Sommerfeld parameter 
space is well described by the boost from the Coulomb 
potential. This part of the parameter space is also the 
part where the range ([TO]) of the potential is less than 
about 0.7 as expected. However, we do not find reso- 
nances in the Coulomb case as we do for the Sommerfeld 
enhancement, because the potential is not localized. To 
learn more about these resonances, we relate the Yukawa 
potential to two model potentials which can be solved 
analytically. 



C. Spherical Well 

The resonance structure of the Sommerfeld- 
enhancement does not appear in the massless limit, 
but we can understand it qualitatively by examining 
the spherical well and relating the well to the Yukawa 
potential [4lj. We look at a potential of the form 
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0, r > L/rrifj,. 
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We set Vo to a value such that the potential integrated 
over volume agrees with the Yukawa potential (pQ). We 
find 



(15) 



where L is the range of the potential in units of ra^ 1 , 
and should be order 1. The Schrodinger equation for this 
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where we are using x = rn^r and 
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This equatio n is solv ed by sines and cosines, so if we 
define p — \/K 2 + e 2 , we can write the general solution 

as 



,x<L 
,x>L 



X~ " " (x) — A sin (px + 7) 



X x> (x) =Bsm(e(x-L) + 5), 



(19a) 
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and using boundary conditions (|5]), we can set 7 = 
and B = 1. We must now match the wave function and 
its derivative at the boundary at x = L. We get the 
equations 



A sin (pL) = sin (5) 
Ap cos (pL) = e cos (6) 

and we solve for the amplitude A. We find: 
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l + (f) 2 cos 2 (pL) 
which results in a boost factor 

i + (f) 2 

l + (f) 2 cos 2 (pL) 

We have plotted this boost factor for L = 1 in figure 
[3] for the same (a, /^-parameter space as before. If we 
compare with figure [TJ we can see the same type of reso- 
nances, although they become suppressed in the Yukawa 
case when the scattering becomes Coulomb-like. The 
structure of ([22]) tells us something about the behavior 
of resonances, even for the Yukawa-case. First off, since 
the denominator is always larger than unity, the maxi- 
mal boost must be of the order ~ K 2 /e 2 . Secondly, if 
the cosine is of order unity, the boost will be order unity 
as well. So the resonances occur when pL/n ~ n + 1/2. 

If we are interested in resonances giving S of order 100 
or more, we can make the approximation e/K <C 1, which 
allows us to use the zeroth order approximation p « K. 
Using the definition of we find 
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so the position of the resonances depends only on the 
ratio of a and /. This suggests that we introduce the 
rotated coordinates 
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(24b) 
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FIG. 3: Boost factor in log 10 for the spherical well with v — 
150km/s. The depth of the well is related to the strength of 
the corresponding Yukawa potential. 



-3.5 1 



o -5.5 1 

II 



© -6.5 1 

-zl 

-7.5 1 
-8 1 



I 



-0.5 



0.5 




log 10 (u)=log 10 (a/f) 



FIG. 4: Boost factor for the spherical well with v = 150km/s 
in the new coordinates u and v. 



Using 



we find that the resonances happens at 
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Since the distance between resonances increases only by 
n 2 , it is clear that they will happen closer and closer 
when viewed in a logarithmic plot, just as we saw in the 
Yukawa-case. The other parameter v controls the order 
of the boost factor, since we have 
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We have shown the boost plot in the rotated coordi- 
nates (?i, v) in figure EJ We suspect that this choice of 
coordinates would also work well in the Yukawa case, so 



2.5 



1.5 



0.5 



-1 



-0.5 



0.5 



1.5 



log 10 (u) = log 10 (o/f) 



2.5 



3.5 



FIG. 5: The Sommerfeld boost factor with v — 150km/s in 
the new coordinates u and v. We have cropped the plot at 
the (a, /) limits from figure [T] 

we calculate the Sommerfeld boost factor in the (u, v) 
coordinates, and the result is shown in figure [5] As is 
clear from the plot, u alone determines the position of 
the peaks, just as for the spherical well. This is in agree- 
ment with [19[, as our coordinate u is identical to their 
parameter e^. 



D. Spherical Slope Well 

We now examine another model potential, the spher- 
ical slope well, to see how the position of resonances 
change compared to the spherical well. We write the 
potential 
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V[r, -\ 0, r> L/ m<j> , 
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and as before, we set Vb to a value such that the potential 
integrated over volume agrees with the Yukawa potential 
(PQ). In this case we find 



Vo = 



(28) 



where L is the range of the potential in units of m^ 1 , 
and should be order 1. The Schrodinger equation for this 
potential is 
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where we are using x = ra^r, K = y jji an d 
2 . If we now do the substitution £ = 
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K 4 / 3 [K 2 x — p 2 ) , we get the Airy equation 
d 2 x 
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which has the general solution 

X = ciAi(0+C2Bi(0- 



(29) 
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We are interested in x a ^ two different values, £o 
£(x = 0) and £l = £(x = L). These are given by 



^ = K~^ (-p 2 ) = -if- 4/3 . 
Cl=^" 4/3 (K 2 L-p 2 ) 



p 2 (31a) 
-if" 4 / 3 e 2 . (31b) 



Since we have £ < in both cases, it makes sense to 
introduce the variable 
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and reexpress the solution x an d its derivative in terms 
of Bessel functions. We have: 



X(0 = [-AJ 1/3 (C) + BJ_ 1/3 (C)] (33a) 
dX, 



d£ 
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where the new coefficients A and B are related to the old 
ones by 

A= - Cl /3 + c 2 /V3 B = Cl /3 + c 2 /V3. (34) 

Since we require x(£o) = 0, equation (|33ap can be used 
to relate A and B, since we have 



A/ 1/3 (Co) = £?J_ 1/3 (Co) => 

B = -p^-T^rrA = tA. 

J-l/3(Co) 



(35) 



At this point we have one free parameter A left to deter- 
mine from the solution x x<L '• The solution outside the 
well is just as before, (|19b|) , and contains one free pa- 
rameter, 5. Both parameters are fixed by matching the 
solution and its derivative at the boundary; we have two 
equations with two unknowns: 



sin(J) =x x<L (x = L) 

dx x<L 
ecos(.5) = — (x = L), 
ax 



(36a) 
(36b) 



and the interesting parameter A is most easily found by 
inserting the right hand sides in the well known trigono- 
metric identity 

e 2 = e 2 sin 2 ((5) + (ecos((5)) 2 

= e 2 (-^ L )A 2 [- J 1/3 (Cl) + ^-1/3 (Cl)] 2 + 
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FIG. 6: Boost factor in log 10 for the spherical slope well with 
v — 150km/s. The depth of the well is related to the strength 
of the corresponding Yukawa potential. 

Using equation (O, we can find the boost factor S: 



S = ^K-^ 3 A 2 [J_ 2/3 (Co) + ^ 2/3 (Co)] 2 

5 _P 4 [J-2/3(Co)+*J 2/3 (Co)] 2 

e4 [-Ji/ 3 (CL)+tJ-i/ 3 (a)] 2 + ' 
+ [j_ 2/3 (a)+tj 2 / 3 (a)] 2 



(38) 



This is a rather complicated expression, but the numer- 
ator can be simplified by using the series expansion for 
J(z) and the following identity for the gamma function: 



r(i/3)r(2/3) = r(i/3)r(i - 1/3) 

= 7rcsc(7r/3) = 2tt/V3 
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Having these, we can prove the neat identity 
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(40) 



Reinserting t, the numerator of equation (|38p can now be 
written as 

num. = [j_ 2/3 (Co) + tJ 2 / 3 (Co)] 2 
= J_ 1/3 (z)- 2 x 

x [J_i/ 3 (Co)J- 2/3 (Co) + Ji/ 3 (Co)J 2/3 (Co)] 2 

(41) 
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Inserting into ([38]) and rearranging the denominator, 
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we find: 
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We have plotted this boost factor in figure [6] in the (u, v)- 
coordinates we introduced for the spherical well. We 
see the same behavior as before, but the expression for 
the boost factor is still too complicated for us to deduce 
where the resonances will be. Let us take another look at 
the denominator of equation (|^2j) . We want to estimate 
the numerical value of the two arguments, Co and Cl- We 
got: 



Co - 3 (-&) - ^ 
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If we are interested in regions of possibly large boost fac- 
tors, we may assume K > e as we did before, and we 
conclude that Co ^ 1 and Cl <C 1. It turns out, that to 
a reasonable approximation, we can set the Bessel func- 
tions equal to their asymptotic values. If we remember 
that 



(g/gr 

I> + 1)' 



Jv{z) — \ — cos (z - vtt/2 - 7T/4) , z > 1 

V 7TZ 

we can identify the two largest terms in the denominator 
of equation (j32j). We are looking for the largest negative 
powers of Cl, because they will give the largest contribu- 
tion to the denominator. We have one term behaving as 

4/3 

~ Cl an d one term from the crossproduct is behaving 
as ~ Cl 1 ' an d the rest of the 6 terms has more positive 
powers. These two worst terms multiply J_!/ 3 (Co), so we 
suspect to have a resonance pattern that follow the zeros 
of J- 1/3 (Co), that is: 
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Thus, we have derived an analytic equation for the po- 
sition of resonances, just as we did for the spherical well, 
and this expression agrees exactly with the resonances in 
figure [6] 



E. Hulthen potential 

Comparing equation ([25]) and ([45]) , we notice that they 
look very similar. Only the fraction in front and the 
'phase' is different. It seems likely that we may fit the 
peaks of the resonances in the Sommerfeld case by an 
expression: 



Lit 2 (n + bf 



(46) 



and this is indeed the case. We find the values L = 0.1592 
and b = 1.006. The value of b is very close to 1, hinting 
that a similar treatment is doable in the Yukawa case, 
yielding a boost factor which depends only on a sine to 
lowest order. A model potential called the Hulthen po- 
tential admits an analytic solution for the s-wave case[42j, 
and the Sommerfeld boost coming from this potential was 
studied in [2l| and later in [22|. 
The potential looks like 



Abe 



-Sr 



1 



-Sr ' 



(47) 



where A and 5 are parameters. We stress that this is 
not a general version of the Yukawa potential, but it is 
a model potential just like the spherical well and the 
spherical slope well. However, unlike these two poten- 
tials, the Hulthen potential reproduces the 1/r-behavior 
of the Yukawa potential in the limit r — »> 0, and it decays 
exponentially instead of having a fixed range. By a pro- 
cedure similar [43] to ours, A and S is found to be A = a 
and S = krricf) = 7r 2 /6. In our notation, the I = case 
of the Sommerfeld boost just before equation 44 in (2lj , 
can be written as 



cosh ^27r^y^ — cos (^2tt 



kf 



>)* 



(48) 



Equation ([^8]) can be derived easily under the assump- 
tion kaf < /3 2 , (which is usually not the case), however, 
a more careful derivation reveals this equation to be true 
in any case. We have plotted the Hulthen boost factor 
in figure and it looks very similar to the Yukawa case 
in figure [U Inspection of equation ([48]) reveals the reso- 
nance pattern to be 



= k(n + l) z 
4(n + D*. 



(49) 



However, for large n, the resonance pattern is shifted 
completely because of the slight difference between 7r 2 /6 
and the fit value. For a velocity of 150km/s we found 
equation ([48]) to be within 10% of the numerical solution 
at 82% of the plotted parameter space and within 30% 
at 98% of the parameter space. For a velocity of lOkm/s, 
the Hulthen boost factor is within 10% in only 43% of the 
parameter space and within 30% in 74% of the param- 
eter space. By using k = 0.15927T 2 and simultaneously 
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FIG. 7: Sommerfeld boost factor in log 10 for the Hulthen 
potential at a relative velocity of 150km/s. 



correcting the phase by 0.006 by hand in equation (jUJ), 
according to the fit parameters, we did somewhat better. 
The modified formula reproduced the resonance pattern 
perfectly, but the agreement of the boost magnitude at 
the first peak became worse. 



III. CONSTRAINING SOMMERFELD 
ENHANCEMENT WITH THE CMB 

There has been some work on constraining models in- 
corporating Sommerfeld Enhancement using the CMB. 
WIMPs annihilating after recombination may contribute 
to reionisation, thereby changing the predicted optical 
depth which can be inferred from the CMB anisotropy 
spectrum. This was studied by [23|. Earlier, during the 
recombination phase, WIMP annihilations with Sommer- 
feld Enhancement also modify the standard scenario. An 
upper bound on this effect can again be inferred from the 
anisotropy spectrum, as was done in (23430| . We will fo- 
cus on WIMP annihilations happening even earlier, in 
the redshift range of approximately 1100<z<2.1-10 6 . 
Annihilations occurring in this redshift range will not in- 
fluence the anisotropy spectrum, but they will distort the 
Black Body spectrum. 

Common to the above mentioned bounds is the fact 
that they use the expected cross section for a thermal 
relic, (crv) TH ~ 3 • 10 _26 cm 3 /s, and multiply this with 
the boost factor to obtain the effective annihilation cross 
section. As was pointed out in |3lj |. this is not strictly 
correct, since the freeze out process is modified by the 
boost mechanism. So to consistently probe the parameter 
space, we must first find the relativist ic annihilation cross 
section ctq, which gives the correct relic abundance. 



A. Relic density calculation 

To catch the full nature of thermal freeze out in models 
with Sommerfeld enhancement, we must do a full calcu- 
lation of the integrated Boltzmann equation because of 
the non-trivial velocity dependence. It should be noted, 
that we do expect something of the same order of mag- 
nitude as (av) TH , since the Sommerfeld boost factor is 
very close to 1 at the time of freeze out. We start from 
the integrated Boltzmann equation for Dark Matter, 

a-*^2^ = (* ann v){nl eq -nl}, (50) 

where a is the scale factor, n x is the number density of 
Dark Matter, n X)eq is the Dark Matter equilibrium num- 
ber density and (cr ann v) is the thermally averaged annihi- 
lation cross section. We normalize the number density to 
the total entropy density s oc a~ 3 by introducing Y = ^f-: 

c\Y 

— = s(a ann v){Y e 2 q -Y 2 } (51) 

We also want to substitute the time parameter t by the 
dimensionless evolution parameter x = -S*-. After some 
manipulations, we find the final form of the Boltzmann 
equation: 

The manipulations leading up to equation ([52]) can be 
found in section [B] of the appendix. In our case, a = 
ctqS(v, . . .), where <Jo is the s-wave annihilation cross 
section which is independent of velocity, and S is the 
Sommerfeld boost factor, which depends on velocity and 
the model parameters. As was noted by the authors 

of [H, [33j| , do is not a completely free parameter but is 

2 

related to a and m x by dimensional analysis, <Jo ~ 

However, the exact factor is highly model dependent, and 
for that reason we keep ao as a parameter. 

Inserting the Sommerfeld cross section in equation ([52]) 
leads to our final equation: 

^ = ^hh (S(v,...)v){Y e l-Y*} (53) 

= \( x ){Y e \-Y 2 }. (54) 

The parameter ao is then found by imposing the bound- 
ary condition that the WIMP must make up all of Dark 
Matter. In this article we are considering the case where 
the WIMP is not its own antiparticle, so we impose the 
boundary condition Udm = 2£l x ,o assuming that there is 
no asymmetry between x an d its antiparticle. 

The equilibrium number density Y eq can be calculated 
exactly under the assumption that Dark Matter follows 
a Maxwell-Boltzmann distribution. We got 

45 r 2 

Y m = gi — — K 2 ( X ), (55) 
47T* g*s 
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where K2 is the modified Bessel function of the second 
kind. The derivation of equation ([55]) can be found in sec- 

tion[B]of the appendix. The values of g*s in equation ([55]) 

1 

and hi in equation ([53]) depends on the temperature, and 
it is found by interpolation in a precomputed table. The 
number of internal degrees of freedom, was set to 2. 
We also need the thermally averaged cross section. If 
we approximate the distribution function for the Dark 
Matter gas by the non-relativistic Maxwell-Boltzmann 
distribution 



(56) 



f(0) = ^xip 2 e~^ 
the thermal average of S becomes 



Jlxi£s(f},...)?e-i*e\ (57) 



This distribution is only normalized to 1 when the upper 
limit goes to infinity, but if we are in the non-relativistic 
limit this is a small correction. (In the calculation we set 
the upper limit dynamically to 4 times the position of 
the peak of the distribution.) Before freeze out the gas 
is not strictly non-relativistic and we should instead use 
the relativistic distribution: 



-x 7 (/3) 



(58) 
(59) 



where 7 is the usual gamma factor and K2 is the modified 
Bessel-function of the second kind. The thermal average 
then becomes 



(Sv) 



k 2 (x) 



5«2 -x( 7 -l) 



(60) 



where fe(^) = K2(x)e x . When 7 becomes close to 1 
we must use the series expansion to calculate 1 — 7 for 
numerical stability. The difference in the required cross 
section ctq from using ([60]) compared to ([57]) was negligi- 
ble, however. 

The numerical solution of equation ([54]) is not entirely 
trivial, especially not since the equation is stiff and we 
need to solve it repeatedly. In section [B] of the appendix, 
we have described the numerical scheme we use for this 
problem, which is the same used in the DarkSUSY soft- 
ware package [34|. We recommend this scheme to others 
interested in Freeze-Out calculations. 



B. Distorting the Black Body spectrum 

When energy is injected into the CMB photons, 
two types of processes are needed to restore a black 
body spectrum: Number changing processes and equi- 
librating processes. Double Compton scattering and 
bremsstrahlung belongs to the first category, while 



Compton scattering and inverse Compton scattering be- 
longs to the second. Double Compton scattering freezes 
out at zdc — 2.1 -10 6 and Compton scattering freezes 
out at zq — 5.4 • 10 4 . Energy which is deposited in the 
photon gas after z^c but before zq will be redistributed 
to give an entropy maximising Planck spectrum with a 
chemical potential fi. Energy input after zc, (but before 
recombination), can not equilibrate and will result in a 
Compton-y distortion [35| of the Planck spectrum. 

The relevant quantity for deriving the size of both ef- 
fects is thus the relative energy input to the CMB during 
these two epochs. Following [36|, we write 



$Pl 

Pi 



t 2 



Pa 



Pi 



-dt 



=1 

Ju 



* 2 2 Fm x (cr ann t>) n\ 

Ay,0(l + ZY ' 



(61) 



where p 1 is the energy density of the CMB photons and z 
is the redshift. F denotes the fraction of energy which is 
transfered to the CMB photons. This is independent of 
redshift when z > 2500, and according to table 1 of [30|, 
this is between 30% and 90% depending on the annihi- 
lation channel. The factor of 2 in equation ([6T]) stems 
from the fact that we, in consistence with our freeze out 
calculation, assume that \ is n °t its own antiparticle. 
We introduce the following relations: 



H 2 {t) 

dt 

dz 



4tt 3 



45m 2 a " "' 
1 45 m 



g.TUl + z) 4 



, 2 >-%-(l + ,)-3 
n x>0 (l + z) 3 = ^{l + zf 



rriy 



87rm v 



pi = P1A 1 + z ) A 

T 7 = T 7> o(l + 20 

Using equations ([62]) in equation ([6T]) yields: 
405 fb 



(62a) 
(62b) 

(62c) 

(62d) 
(62e) 
(62f) 



bpi 
Pi 



64tt 5 

z(ti) 

(*2) 

405 



9, 2 F 
7r m 



^avg(^ 7 • • •) 



7,0 



1 



-dz 



647T 5 V 7T 



l + z 
lOOGeV 



ctq 



10- 26 cm 3 /s 



x (n x hf c_ 7 / s avs 

Jz(t 2 ) 



(*,...) 



l + z 



-dz. 



(63) 



Here C-j ~ 2.8696 • 10 _T is a numerical constant. Since 
the Dark Matter is very cold at this time, the distribu- 
tion function is strongly peaked around its mean value, 
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FIG. 8: Required annihilation cross section ao in units of 
10 -26 cm 3 /s to explain the total Dark Matter abundance for 
a 200GeV WIMP and a kinetic decoupling temperature of 
8MeV. 

so to a good approximation, we may write S avg (z) 
S (/3 meSin (z)) . We may also assume that the Sommerfeld 
enhancement has saturated at this stage, making the ap- 
proximation S(f3 mean (z)) ~ #(/3 me anO0l))), allowing US 

to pull S outside the integral which can then be done 
analytically: 



Pi 



405 
64tt 5 



9**F 



lOOGeV 



(O x /i) 2 C_ 7 5(/3 mean (zi))ln 



10-26 cm 3/ s 



1 + z(t 2 ) 



(64) 



In these calculations we have assumed a radiation dom- 
inated universe as described by the Friedmann equation 
(|62ap , which breaks down at z ~ z eq ~ 3300, well before 
recombination. But as can be seen from the analytical 
approximation, the dependence on zfo) is logarithmic, 
so the error in doing this is insignificant. When plotting 
the \i and \y\ -distortions, we have put F = 1 for conve- 
nience, but since the dependence is linear in F, it can be 
reinstated by multiplying the distortions by F. 



C. Results 

For each point in the (a, /) -parameter space, we solved 
the Boltzmann equation from x\ — 1 to x 2 = m x /T 7) o to 
find the value of <r . This depends weakly on the WIMP 
mass as well as the kinetic decoupling temperature but 
the overall dependence on the Sommerfeld parameters 
can be seen in figure [8] The kinetic decoupling tempera- 
ture was taken as a parameter, and the effect on the freeze 
out process is shown on figure [9] We find that this effect 
can be as large as 30% in agreement with [3l|, [37j . We 
expect ao to be nearly independent of the WIMP mass 
and this is confirmed by figure [10] which shows maxi- 
mally one percent difference between a 200GeV WIMP 
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FIG. 9: The ratio a KD=8MeV /(j KD=500MeV for a 200GeV 
WIMP. 
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FIG. 10: The ratio (tq/gq of cross sections for a 200GeV 
WIMP and a lOOOGeV WIMP respectively, which decouples 
kinetically at xkd = 2 • 10 3 . 

and a lOOOGeV WIMP for a kinetic decoupling value of 
xkd = 2 • 10 3 for both. However, in a real particle physics 
model we expect xkd to have a slight mass dependence, 
so the actual effect of different masses may be somewhat 
bigger [38]. 

For each point in the parameter space, we solve the 
integral in equation ([63]) on both the interval which is 
relevant for Compton-|?/| distortions as well as the inter- 
val relevant for the CMB photons to develop a chemical 
potential /i. The results for a 200GeV particle are shown 
in figure [11] and [12] Considering the current bounds on 
li- and ^/-distortions from FIRAS Q of |/i| < 9-10" 5 , 
| y | < 1.5-10 -5 , only a small portion of the parameter 
space can be ruled out. But there is a rather large part 
of lower right part of the parameter space which is close 
at saturating the current bound. The analysis for the ji- 
distortion was already carried out in [37| , and our results 
agrees with theirs. 

As is evident from figure [11] and [121 the two bounds are 
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FIG. 11: The magnitude of the Compton-|?/| parameter in 
log 10 for a 200GeV WIMP and a kinetic decoupling temper- 
ature of 8MeV. 
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FIG. 12: The magnitude of the chemical potential \i in log 10 
for a 200GeV WIMP and a kinetic decoupling temperature of 
8MeV. 

degenerate in the way that they both tends to rule out 
the resonances and the lower right part of the parameter 
space. With the current limits, the \y\ -bound is always as 
strong or stronger than the /i-bound. One would suspect 
this to be the case, since the ^-distortion happens at a 
later time than the /i-distortion, giving the WIMPs more 
time to cool. This is indeed true, but only for a small 
subset of the parameter space, at the resonances in the 
lower left part. For the rest of the parameter space, the 
Sommerfeld enhancement has already saturated at this 
point, and no further enhancement is possible. 

We can consider what would be the allowed possibil- 
ities for the annihilation cross section in a halo having 
a fiducial velocity dispersion of 150km/s by removing all 
points that exceeds either bound. We have plotted this in 
figure [13] for our 200GeV example WIMP. It is suspected 
that a new FIRAS-like satellite, if built, could bring the 
bound on | y | down to the order ~ 10 -7 . In figure [TH we 
have shown what figure [13] would look like with a fidu- 
cial bound of \y\ < 10 _T . It is also worth comparing this 
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FIG. 13: Logarithm of the annihilation cross section in units 
of 10 _26 cm 3 /s in a halo with velocity dispersion of approxi- 
mately 150km/s. White points are ruled out by either fi- or 
|?/|-bound. 




log 10 (oc) 

FIG. 14: Logarithm of the annihilation cross section in units 
of 10 _26 cm 3 /s in a halo with velocity dispersion of approxi- 
mately 150km/s. White points can be ruled out if the bound 
gets improved by a factor of 500. 



bound with the anisotropy bound. From [3C 
the bound 



we have 



(cTann^) 



saturated 



< 



360 • 10- 26 cm 3 /s 



ITeV' 



(65) 



where F YC is the average fraction of energy being trans- 
fered to the CMB at the time of recombination. This has 
different values depending on the annihilation channel, 
as can be seen in table 1 of [30j |. However, it is roughly 
of order 30% for most processes. 

As we discussed earlier, plugging in the standard value 
(^ann^) T H f° r a thermal relic is not strictly accurate. It 
is easily fixed however, by using our calculated values of 
ctq, as well as the thermally averaged Sommerfeld factor 
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FIG. 15: Logarithm of the annihilation cross section in units 
of 10 -26 cm 3 /s in a halo with velocity dispersion of approx- 
imately 150km/s. White points are ruled out by the CMB 
anisotropy bound. 



at recombination. We get 



10- 26 cm 3 /s 



-^avg 



. . 360 m x 

{Zrc '"' }< FlTTeV' 



(66) 



In figure [15] we have again showed the effective boost fac- 
tor, but this time with the CMB bound, equation ([66]) . 
It is clear, that this bound is stronger than even the fore- 
casted 7/-bound by roughly an order of magnitude. 



IV. CONCLUSION 

We have analysed the Sommerfeld enhancement mech- 
anism in detail and found coordinates in which the po- 
sition of resonances depends only on one coordinate, u, 
and we have derived analytical expressions for the posi- 
tion of resonances for two model potentials. A numeric 
treatment of the Sommerfeld case showed that the a sim- 
ilar simple relationship exists for this potential, and we 
found the following fit to agree excellently: 



0.1592 (n + 1.006) 2 



0,1, 



(67) 



This is a nice result, since knowing the position of reso- 
nances beforehand is helpful for doing numerical calcula- 
tions. 

In the second part we did the full freeze out calculation 
for the Sommerfeld parameter space by solving the Boltz- 
mann equation all the way through chemical and kinetic 
freeze out, and we then calculated the effect of annihila- 
tions on the CMB blackbody spectrum. We found that 
only a small portion of the parameter space, directly on 
the resonances, can be ruled out by the current bound 
from FIRAS. We also noted that a future measurement 
of the \y\ — or fi— distortion would rule out a huge portion 
of the parameter space. But as we showed, the anisotropy 



bound is already stronger than this and is likely to im- 
prove with Planck data. However, we think that it is 
worth mentioning, that the \y\— and fi— bounds are ob- 
tained at a different epoch and by different observations 
than the anisotropy bounds, and thus should be consid- 
ered complementary to those. 



Appendix A: Sommerfeld enhancement for the 
Coulomb potential 

We want to derive equation (fT3]h the boost factor for 
the Coulomb potential. We start from the radial equation 
(j5j) with = and do the substitution x = am x r: 



j_fx 

m x dr 2 

fx 
dx 2 



/ a 
\ r 



rn 



tP 2 ). 



4 



(Al) 



x 



00. 



(A2a) 
(A2b) 



We analyze equation (|A1|) the usual way by considering 
the asymptotic limits of the equation. We got: 

fx 

dx 2 

fx 

dx 2 

The analytic solution of (|A2b|) is just an exponential, 
while the analytic solution of (|A2b|) is more complicated. 
In general we have 

X (X) = Ciyf{x)Ji{2y/x) + C 2 y/xYi{2y/x), 

where J\ and Y\ are the Bessel functions of the first and 
second order, respectively. However, y/xYi(2y/x) is not 
well behaved for x — » 0, so we are left with the C\ term. 
Since we are looking at x <C 1, we may use J\{z) ~ z/2. 
Our Ansatz for the solution then becomes: 



X 



xe 



c v(x), 



(A3) 



where v(x) interpolates between the asymptotic solu- 
tions. Taking the second derivative of the Ansatz (|A3|) 
yields 

= e' ecx [{2ie c - e 2 c )v + (2 + 2ie c x)v' + xv"] , 

which can be inserted into equation (jAip . The prime on 
v denotes differentiation w.r.t x. The result is 







+ (2 + 2ie c x)v' + (2ie c + l)v. 



(A4) 



We can bring this equation on a more recognizable form 
by making the substitution z = —2iecx. This leads to 







fv 

''dz 2 



+ (2 - z) 



dv 
~dz 



1 



2ec) 



(A5) 
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which is the confluent hyper geometric equation. The 
solution which is regular in x = z = is then: 



(A6) 



v(z)=CF[l-—,2,z). 



We now need to apply the asymptotic boundary condi- 
tion (|8bp to our solution, so we need the asymptotic for- 
mula for F(a,b,z) in the limit x — » oo or, equivalent ly, 
z — > — zoo. To do this, define the following: 



g(a,b,z) = 



(c) n = c(c + l)(c + 2) • • • (c + n - 1), (c)o = 1 



F(a,M) = 



r(6) 



r(6-a 

r(6) 



-#(a, l + a-6, -z)(-z)" a + 



r(a) 



#(6 - a, 1 - a,z)e 2 z a b . 



Since g(a, 6, z) — > 1 as |z| — > oo, F(a, 6, z) has the asymp- 
totic form 



F°°(a,M) 



r(6) -hri^. 



T(6 - a) 



r(a) 



Now consider the gamma- function T(b — a) for our values 
of a = 1 - i/Oc) and 6 = 2: 

r (b - a) = r (2 - (1 - ) = r (a*) = r (a)* , 

and if we now define 77 by T (a) = |T (a)| e 2?7 , we find the 
following: 



r(6) r(6) 



T(6-a) |r(a)| 

r(6) _ r(6) 



I» |I»| 
The asymptotic form F 00 can now be written like: 

r(6) 



(A8a) 
(A8b) 



F°°(a,b,z) 



ir n(_ y \-a _±_ e -ir]+z z a-b 



|r(a)| 



(e"»(-z) _0 + 



x sin ( ecx H ln(2ec#) + rj ) • 

Substituting the asymptotic solution back into the 
Ansatz (|A3j) yields 



x°°(a:) = Ci- 



- (l-5fe)^ 

x sin (e c x + ^— \n(2e c x) + rj ) , (A9) 



where Ci is now set by the boundary condition (j8b| . We 
find 



Ci = 



(A10) 



which we can finally insert into (jA6|) and the Ansatz (jA3|) : 



xF(l-— ,2,-2ie c x). (All) 



We are interested in the boost factor (0), so we calculate 





2 


k 





^dx 

ec dx 



(x = 0) 



(A12) 



We take the derivative of our solution (jAll|) and evaluate 
it in x = 0: 



dx 



{x = 0) = e c r 1 



1 - 1 e _ir ^ 



2e c 



sinh (-^) 
1 - e- 71 "/^ ' 



e ze c 



(A13) 



which is the result obtained by Sommerfeld and quoted 
in 0. 



Appendix B: How to solve the integrated Boltzmann 
equation 

We want to derive the evolution equation ([52]) for the 
number density of Dark Matter. Starting from the inte- 
grated Boltzmann equation, we get: 

a _ 3 d{ngZ) =<(Tannt; ){ n 2 jeq _ n 2} > (B1) 

where a is the scale factor, n x is the number density of 
Dark Matter, n x , e q is the Dark Matter equilibrium num- 
ber density and (cr ann v) is the thermally averaged annihi- 
lation cross section. We normalize the number density to 
the total entropy density s oc a~ 3 by introducing Y = — : 



dr 



s(a ann v){Y e \-Y 2 }. 



(B2) 
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We want to substitute the time parameter t by the di- 
mensionless evolution parameter x = ^ . We need a few 

formulas to proceed. Using s oc a -3 , we find: 



dsa 3 
~dT 



= ^ s = -3sH. 



The entropy density s is given by 

s(T) = kg* s {T)T\ 



(B3) 



(B4) 
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where /c = is a constant and g*s is the relativistic 



degrees of freedom for the entropy density. Taking the 
time derivative yields: 



ds • 
s = —T = T 
dT 



' dT 



T + 3kg* s T 



= 3s — 



dg.. 



3l 

T 

i i 



35 



T gl hi 
T g*s 



(B5) 



where the last line defines the parameter hi . Equating 
equation (|B3|) and (|B5|) , yields 



T 
T 



-H 



g*s 
i i 
g^. ft* 



1 



TYlph^; 



(B6) 



where the Friedmann equation in the radiation domi- 
nated universe has been used in the last equation. Using 
the chain rule on the left hand side of (f5T|) with this 
equation gives 



dY _ dY dx dT _ dY T 
~dt ~dxdT~dt ~dx X \ T 



dY Utt s /m x \2 1 

= WlHv) ^ (T) ^T- (B7) 

This gives our final form of the evolution equation: 
1. Calculating the average number density 



The number density of Dark Matter in equilibrium, 
n eq , can be found by integrating the distribution function 
/(x, p) over momentum space. We assume a Maxwell- 
Boltzmann distribution with zero chemical potential: 



/(x,p)=^e t . 



(B9) 



Using E(p) = \/p 2 + m 2 , and introducing x = tu x /T 
and a = p/T we get: 



/ 



d 3 P 



e T 



4ngT 3 f 



/OO 



(2tt) 2 
(2tt) 3 



271"^ £ 



47T 4 



(BIO) 



where K2 is the modified Bessel function of the second 
kind. 



2. Solving the differential equation 

Equation ([54]) is an example of a stiff differential equa- 
tion, meaning that it involves vastly different scales. In 
this case it is roughly the timescale at which annihila- 
tions take place versus the timescale at which expansion 
happens. Explicit solvers like Runge-Kutta either fails 
or need a very low step size to maintain stability. We 
implement the implicit scheme, which is also used in the 
software package DarkSUSY 0, Eo|. The derivative y' 
of a function can be approximated as follows: 



Forward difference 



(Blla) 



y' i+1 ~ — — Backward difference. (Bllb) 



Equation (|Bllb|) inserted in the general ODE y' = 
f(x,y) gives rise to the implicit backward euler scheme: 

Vi+i =Vi + hf i+1 (x, y) + (ft) . (B12) 

Adding equations ([Blip yields the implicit trapezoidal 
rule 



y i+1 =Vi + \ (fi+i +fi) + G {h 2 



(B13) 



These two schemes are the s = and s = 1 Adams- 
Moulton methods respectively. We can use the difference 
between the two methods as an error-estimate for con- 
trolling the step size. An implicit scheme like this usu- 
ally results in a system of non-linear algebraic equations 
which must then be solved numerically in each step. For- 
tunately, for the Boltzmann-equation, and using method 
(|B12|) or (|B13|) , this is just a second order equation which 
can be solved analytically. We introduce the variables 
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suggested in [34[: 



Y eq q 


(B14a) 


u h\ i+ i 


(B14b) 


K 


(B14c) 


c = 2Y l + u[(qt +1 +pqt)-pY?] 


(B14d) 




(B14e) 



The Euler method (|B12j) for the Boltzmann equation 



is 



If we let denote the Euler estimate of the next step, 
we estimate the relative error err as 



err = 



Yi+l 



(BIT) 



4Y i+1 =4y i ^4h\ i+1 (ql 1 -Y^ 1 ) 



Auy- 



i+l 



Yi+i — 



Yi 



1 ip y/l + UC! 

2u 

d 



+1 



1 



21 + vT" 



(B15) 



where we have chosen the solution which gives a positive 
Y i+ i. It works the same way for the trapezoidal rule 
SI: 



2Y i+1 = 2Y i + h [Ai(g? - if) + A i+1 (^ +1 - 
= 2Yi + u [pq* + q 2 i+1 - pY 2 } - uY* 



r2 



c-uY^ 



Yi+i 
Yi+i 



-1 ± \/l + ucu 



1 + Vl + wc 



(B16) 



and since the trapezoidal rule is second order in h, we 
modify the step size according to 



hnext = min ( hS 




5h 



(B18) 



where S is a safety factor set at 0.9 and eps is the 
wanted accuracy. We have also demanded that h can 
only grow with a factor of 5 in each step. 
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